rm(list = ls())
library(DropletUtils)
## Warning: package 'DropletUtils' was built under R version 4.1.2
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## Warning: package 'matrixStats' was built under R version 4.1.2
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## Loading required package: GenomicRanges
## Warning: package 'GenomicRanges' was built under R version 4.1.2
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Warning: package 'GenomeInfoDb' was built under R version 4.1.2
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
library(scater)
## Loading required package: scuttle
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.1.2
library(scran)
## Warning: package 'scran' was built under R version 4.1.2
library(ggplot2)
library(Seurat)
## Warning: package 'Seurat' was built under R version 4.1.2
## Attaching SeuratObject
##
## Attaching package: 'Seurat'
## The following object is masked from 'package:SummarizedExperiment':
##
## Assays
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.2
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:Biobase':
##
## combine
## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
##
## intersect
## The following objects are masked from 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union
## The following objects are masked from 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union
## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union
## The following object is masked from 'package:matrixStats':
##
## count
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(viridis)
## Warning: package 'viridis' was built under R version 4.1.2
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.1.2
library(EnhancedVolcano)
## Loading required package: ggrepel
## Warning: package 'ggrepel' was built under R version 4.1.2
## Registered S3 methods overwritten by 'ggalt':
## method from
## grid.draw.absoluteGrob ggplot2
## grobHeight.absoluteGrob ggplot2
## grobWidth.absoluteGrob ggplot2
## grobX.absoluteGrob ggplot2
## grobY.absoluteGrob ggplot2
library(formattable)
##
## Attaching package: 'formattable'
## The following object is masked from 'package:scater':
##
## normalize
## The following object is masked from 'package:BiocGenerics':
##
## normalize
library(pheatmap)
library(ape)
## Warning: package 'ape' was built under R version 4.1.2
##
## Attaching package: 'ape'
## The following object is masked from 'package:dplyr':
##
## where
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
## in the original pheatmap() are identically supported in the new function. You
## can still use the original function by explicitly calling pheatmap::pheatmap().
##
## Attaching package: 'ComplexHeatmap'
## The following object is masked from 'package:pheatmap':
##
## pheatmap
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.1.2
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:formattable':
##
## area
library(scales)
## Warning: package 'scales' was built under R version 4.1.2
##
## Attaching package: 'scales'
## The following objects are masked from 'package:formattable':
##
## comma, percent, scientific
## The following object is masked from 'package:viridis':
##
## viridis_pal
library(wesanderson)
library(ggrepel)
library(ggeasy)
## Warning: package 'ggeasy' was built under R version 4.1.2
library(data.table)
## Warning: package 'data.table' was built under R version 4.1.2
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following object is masked from 'package:SummarizedExperiment':
##
## shift
## The following object is masked from 'package:GenomicRanges':
##
## shift
## The following object is masked from 'package:IRanges':
##
## shift
## The following objects are masked from 'package:S4Vectors':
##
## first, second
library(scuttle)
library(SingleR)
## Warning: package 'SingleR' was built under R version 4.1.2
library("Seurat")
library("monocle3")
##
## Attaching package: 'monocle3'
## The following objects are masked from 'package:Biobase':
##
## exprs, fData, fData<-, pData, pData<-
library("org.Hs.eg.db")
## Loading required package: AnnotationDbi
## Warning: package 'AnnotationDbi' was built under R version 4.1.2
##
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
##
## select
##
library("dplyr")
addTaskCallback(function(...) {set.seed(100);TRUE})
## 2
## 2
#Load object
combined_sce1 <- readRDS("FinalObj.RDS")
combined_sce2 <- combined_sce1
sce2 <- as.Seurat(combined_sce2)
Idents(sce2) <- sce2$clusters
#
sce1 <- sce2
sce1
## An object of class Seurat
## 31719 features across 43436 samples within 1 assay
## Active assay: originalexp (31719 features, 0 variable features)
## 2 dimensional reductions calculated: PCA, UMAP
#UMAPs For Figure 5A and 5B
color_samp <- wes_palette("Darjeeling1",3)
wt_col <- color_samp[2]
ko_col <- color_samp[3]
my_color_palette <- hue_pal()(length(levels(sce1$clusters)))
DimPlot(sce1) + ggtitle("UMAP by Clusters")
DimPlot(sce1, group.by = "condition", cols = c(wt_col, ko_col), label = FALSE) + ggtitle("UMAP by WT vs KO") + theme(legend.text=element_text(size=13))
#MAKE HEAT MAP, CLUSTERPLOTS, VLN
markers <- FindAllMarkers(sce1, min.pct = 0.15, min.diff.pct = 0.04, logfc.threshold = .25, test.use = "MAST") #Try with #MAST
## Calculating cluster 1
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 2
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 3
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 4
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 5
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 6
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 7
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
top20 <- markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
sce1 <- ScaleData(object = sce1, block.size=100)
## Centering and scaling data matrix
temp3 <- GetAssayData(sce1, slot = "scale.data")
mat<- temp3[unique(c("JUN", top20$gene)), ] %>% as.matrix()
heatmapanno <- HeatmapAnnotation(Group = sce1$condition, Cluster = Idents(sce1), col = list(Group=c("WT" = wt_col, "KO" = ko_col), Cluster =c("1" = "#F8766D", "2" = "#C49A00", "3" = "#53B400", "4" = "#00C094", "5" = "#00B6EB", "6" = "#A58AFF", "7" ="#FB61D7")))
col_fun = circlize::colorRamp2(c(-1, 0, 2), c("#140b34" , "#84206b", "#f6d746"))
##Figure 5D
Heatmap(mat, name = "Scaled Expression",
column_split = Idents(sce1),
#column_title = "Cluster",
cluster_columns = FALSE,
show_column_dend = FALSE,
cluster_column_slices = FALSE,
column_title_gp = gpar(fontsize = 7),
column_gap = unit(.6, "mm"),
cluster_rows = TRUE,
show_row_dend = TRUE,
row_dend_width = unit(20, "mm"),
col = col_fun,
row_names_gp = gpar(fontsize = 8),
column_title_rot = 0,
top_annotation = heatmapanno,
show_column_names = FALSE,
use_raster = TRUE,
raster_quality = 4)
#Figure 5F sce3 <- subset(sce1, ident = 7, invert = TRUE) v1 <- VlnPlot(sce3, “BAX”, split.by = “condition”, split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0) v2 <- VlnPlot(sce3, “TNFRSF12A”, split.by = “condition”, split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0) v3 <- VlnPlot(sce3, “JUN”, split.by = “condition”, split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0) v1+v2 + v3
##Figure 5E
Idents(sce1) <-sce1$condition
markers1 <- FindAllMarkers(sce1, min.pct = 0.05, min.diff.pct = 0.02, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster KO
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
top <- markers1 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
top$avg_log2FC < -.10, '#F2AD00',
ifelse(top$avg_log2FC > .10, '#00A08A',
'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
apotosis1 <-c("FAS",
"TNFRSF10A",
"TNFRSF10B",
"TNFRSF10C",
"TNFRSF10D",
"TNFRSF11B",
"TNFSF10",
"TNFRSF1A",
"TNFRSF12A",
"FADD",
"CFLAR",
"CASP1",
"CASP2",
"CASP3",
"CASP4",
"CASP6",
"CASP7",
"CASP8",
"CASP9",
"CASP10",
"CASP14",
"NAIP",
"BIRC2",
"BIRC3",
"XIAP",
"BIRC5",
"BIRC6",
"BIRC7",
"BCL2",
"MCL1",
"BCL2L1",
"BCL2L2",
"BCL2A1",
"BCL2L10",
"BAX",
"BAK1",
"BOK",
"BID",
"BCL2L11",
"BMF",
"BAD",
"BIK",
"HRK",
"PMAIP1",
"BNIP3",
"BNIP3L",
"BCL2L14",
"BBC3",
"BCL2L12",
"BCL2L13",
"APAF1",
"CYCS",
"DIABLO",
"HTRA2",
"AIFM1",
"ENDOG",
"CARD8",
"CARD6",
"NOX5",
"TP53",
"CDKN1A",
"CDKN1B",
"CDKN2A",
"CDKN2B")
ofinterest <- c("MKI67",
"TOP2A",
"NES",
"MAP2",
"TH",
"EN1",
"NR4A2",
"HIF1A",
"TP53",
"RELA",
"NFKB1",
"NFKB2",
"JUN",
"FOS", "LMO3", "HES5")
geneofinterest <- c(ofinterest, apotosis1)
samplabels1 <- c("PEG10", "DYNC1I2", "JUN",
"NMT1",
"PPDPF",
"PHPT1",
"NDUFA13",
"NDUFC2",
"NDUFA11",
"NME2",
"NDUFA12",
"RPS27L",
"INSM1",
"FABP5",
"HSPA1A",
"WLS",
"MICOS13",
"ASCL1",
"SOX2")
sce4 <- subset(sce1, ident == 3 | ident == 5 | ident == 6)
Idents(sce4) <-sce4$condition
markers2 <- FindAllMarkers(sce4, min.pct = 0.01, min.diff.pct = 0.01, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster KO
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
top2 <- markers2 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
top2$avg_log2FC > .10 & top2$p_val_adj < .01, '#00A08A',
ifelse(top2$avg_log2FC < -.10 & top2$p_val_adj < .01, '#F2AD00',
'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
samplabels<- c("RPS27L", "PHLDA3", "TUBB2A","TUBB2B","RPL27A","PHPT1","CDKN1A","BAX",
"NMT1",
"BBC3",
"CDKN2B",
"TNFRSF10B")
EnhancedVolcano(top2,
lab = top2$gene,
x = 'avg_log2FC',
y = 'p_val_adj',
title = 'WT vs KO Top DE Genes',
pCutoff = .01,
FCcutoff = .10,
pointSize = 1.0,
labSize = 3.3,
selectLab = samplabels,
col=c('blue', 'green', 'black', 'red3', "yellow", "brown", "grey"),
colAlpha = .8,
colCustom = keyvals,
drawConnectors = TRUE,
widthConnectors = 0.1,
lengthConnectors = unit(0.01, "npc"),
cutoffLineType = "dotted",
cutoffLineCol = "#00000066",
#titleLabSize = 10,
#subtitleLabSize = 8,
#captionLabSize = 8
) + xlim(-.75, .75)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
##Barplot and Gene plot 5SE 5SD 5SB
table2 <- as.matrix(table(sce1@meta.data$condition, sce1@meta.data$cluster))
table2[1,] <- table2[1,] * table(sce1$condition)[2]/table(sce1$condition)[1]
table3 <- data.table(t(t(table2) / colSums(table2)))
colnames(table3) <- c("Condition", "Cluster", "Fraction")
ggplot(data=table3, aes(x=Cluster, y=Fraction, fill=Condition)) +
geom_bar(stat="identity") + scale_fill_manual(values = c(ko_col, wt_col)) + ggtitle("Clusters by Condition")
sce1 <- sce2
Idents(sce1) <- sce1$condition
sce1 <- subset(x = sce1, idents = "KO")
th <- length(WhichCells(sce1, slot = 'counts', expression = TH > 0)) / length(WhichCells(sce1, slot = 'counts'))
mki67 <- length(WhichCells(sce1, slot = 'counts', expression = MKI67 > 0)) / length(WhichCells(sce1, slot = 'counts'))
map2 <- length(WhichCells(sce1, slot = 'counts', expression = MAP2 > 0)) / length(WhichCells(sce1, slot = 'counts'))
all_ko <- c(map2, th, mki67)
sce1 <- sce2
Idents(sce1) <- sce1$condition
sce1 <- subset(x = sce1, idents = "WT")
th <- length(WhichCells(sce1, slot = 'counts', expression = TH > 0)) / length(WhichCells(sce1, slot = 'counts'))
mki67 <- length(WhichCells(sce1, slot = 'counts', expression = MKI67 > 0)) / length(WhichCells(sce1, slot = 'counts'))
map2 <- length(WhichCells(sce1, slot = 'counts', expression = MAP2 > 0)) / length(WhichCells(sce1, slot = 'counts'))
all_wt <- c(map2, th, mki67)
temp <- data.table(
Gene = c("MAP2", "TH", "MKI67"),
Group =c("WT", "WT", "WT", "KO", "KO", "KO"),
Values = c(all_wt, all_ko)
)
ggplot(temp, aes(x = Gene, y = Values, fill = Group)) +
geom_bar(stat = "identity", position = "dodge") + ggtitle("Fraction of Cells Expressing Each Gene") + ggeasy::easy_center_title() +
theme(plot.title = element_text(face="bold")) +theme(text=element_text(size=8)) + scale_fill_manual(values = c(ko_col, wt_col))
##Supplement 5C
sce1 <- sce2
sce1.list = list(sce1[, sce1@meta.data$condition == "WT"], sce1[, sce1@meta.data$condition == "KO"])
f1 <-FeaturePlot(sce1.list[[1]], feature="MAP2", order=T, pt.size=1.0) + ggplot2::labs(title="MAP2 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f2 <-FeaturePlot(sce1.list[[2]], feature="MAP2", order=T, pt.size=1.0) + labs(title="MAP2 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#PBX1
f3 <-FeaturePlot(sce1.list[[1]], feature="PBX1", order=T, pt.size=1.0) + labs(title="PBX1 (WT)") + scale_colour_gradientn(colours =terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f4 <-FeaturePlot(sce1.list[[2]], feature="PBX1", order=T, pt.size=1.0) + labs(title="PBX1 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#MKI67
f5 <-FeaturePlot(sce1.list[[1]], feature="MKI67", order=T, pt.size=1.0) + labs(title="MKI67 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f6 <-FeaturePlot(sce1.list[[2]], feature="MKI67", order=T, pt.size=1.0) + labs(title="MKI67 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#SERT
f7 <-FeaturePlot(sce1.list[[1]], feature="SLC6A4", order=T, pt.size=1.0) + labs(title="SLC6A4 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f8 <-FeaturePlot(sce1.list[[2]], feature="SLC6A4", order=T, pt.size=1.0) + labs(title="SLC6A4 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#TPH1
f9 <-FeaturePlot(sce1.list[[1]], feature="TPH1", order=T, pt.size=1.0) + labs(title="TPH1 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f10 <-FeaturePlot(sce1.list[[2]], feature="TPH1", order=T, pt.size=1.0) + labs(title="TPH1 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
(f1/f2) | (f3/f4) | (f5/f6)
##FIGURE 5C
#load cellline data
cellline <- readRDS("cellline_clean.rds")
cellline <- as.Seurat(cellline)
## Warning: The column names of the counts matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: The column names of the data matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from PC to PC_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to PC_
#La Manno dataset used a reference
combined_sce1 <- as.SingleCellExperiment(sce1)
combined_sce1$clusters <- sce1$clusters
combined_sce_lamanno <- scRNAseq::LaMannoBrainData("human-embryo")
## snapshotDate(): 2021-10-19
## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
## loading from cache
## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
## loading from cache
# SingleR() expects reference datasets to be normalized and log-transformed.
combined_sce_lamanno <- logNormCounts(combined_sce_lamanno)
predictions2 <- SingleR(test=combined_sce1, assay.type.test=2,
ref=combined_sce_lamanno, labels =combined_sce_lamanno$Cell_type, de.method="wilcox", clusters=combined_sce1$clusters)
predictions2$pruned.labels #Didn't assign
## [1] "hNbM" "hNbM" "hProgFPL" "hNbM" "hProgFPL" "hProgFPL" "hPeric"
sce5 <- as.Seurat(combined_sce1)
Idents(sce5) <- sce5$clusters
sce5 <- RenameIdents(object = sce5, `1` = "hNbM", `2` = "hNbM", `3` = "hProgFPL", '4'= "hNbM", "5" = "hProgFPL", "6" = "hProgFPL", "7" = "hPeric")
DimPlot(sce5) + ggtitle("UMAP by Annotation")
combined_sce_lamanno1 <- subset(combined_sce_lamanno , ,Cell_type == c('hNbM') | Cell_type == c('hNbML1') | Cell_type == c('hProgFPL'))
predictions3 <- SingleR(test=combined_sce1, assay.type.test=2,
ref=combined_sce_lamanno1, labels =combined_sce_lamanno1$Cell_type, de.method="wilcox")
combined_sce1$pruned.labels <- predictions3$pruned.labels
sce6 <- as.Seurat(combined_sce1)
v1 <- VlnPlot(sce6, "BAX", group.by = "pruned.labels", cols = c(wt_col, ko_col), pt.size = 0, split.by = "condition")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
v2 <- VlnPlot(sce6, "TNFRSF12A", group.by = "pruned.labels", cols = c(wt_col, ko_col), pt.size = 0, split.by = "condition")
v3 <- VlnPlot(sce6, "JUN", group.by = "pruned.labels", cols = c(wt_col, ko_col), pt.size = 0, split.by = "condition")
v1 + v2 + v3
cellline <- readRDS("cellline_clean.rds")
cellline <- as.Seurat(cellline)
## Warning: The column names of the counts matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: The column names of the data matrix is NULL. Setting cell names to
## cell_columnidx (e.g 'cell_1').
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from PC to PC_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to PC_
f1 <- FeaturePlot(sce1, feature="HES5", order=T, pt.size=1.0) + labs(title="HES5") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
f1 <- FeaturePlot(sce1, feature="HES5", order=T, pt.size=1.0) + labs(title="HES5") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
#ALTERNATE FIGURE REQUESTS
sce4 <- subset(sce1, ident == 1 |ident == 2 | ident == 4)
Idents(sce4) <-sce4$condition
markers2 <- FindAllMarkers(sce4, min.pct = 0.01, min.diff.pct = 0.01, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster KO
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
top2 <- markers2 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
top2$avg_log2FC >= .10 & top2$p_val_adj < .01, '#00A08A',
ifelse(top2$avg_log2FC <= -.10 & top2$p_val_adj < .01, '#F2AD00',
'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
samplabels<- c("PEG10",
"NMT1",
"HSPA8", "RPS27L",
"NDUFA12",
"DYNLT1",
"PDCD5",
"PHPT1",
"KRTCAP2", "SYT14",
"SOX2",
"HSPA1B",
"HSPA1A","HSPA2A","DNAJB1","HSPA8")
EnhancedVolcano(top2,
lab = top2$gene,
x = 'avg_log2FC',
y = 'p_val_adj',
title = 'WT vs KO Top DE Genes',
pCutoff = .01,
FCcutoff = .10,
pointSize = 1.0,
labSize = 3.3,
selectLab = samplabels,
col=c('blue', 'green', 'black', 'red3', "yellow", "brown", "grey"),
colAlpha = .8,
colCustom = keyvals,
drawConnectors = TRUE,
widthConnectors = 0.1,
lengthConnectors = unit(0.01, "npc"),
cutoffLineType = "dotted",
cutoffLineCol = "#00000066",
#titleLabSize = 10,
#subtitleLabSize = 8,
#captionLabSize = 8
) + xlim(-.8, .8)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.